function [ v, zeta , Lh , L , vL , C , ce , ...
           pop , f , w , wh , u , H , h , s ] = SolveEq( Grids , Params , GE )

% solve model       
[ v , zeta , Lh , L , vL , C ] = ...
    SolveC( Grids , Params , GE ) ;

% Recover actual population density
L = Lh / sum( Lh .* Grids.mx ) ;
pop = L .* Grids.mx ; pop = pop / sum(pop) ;

% construct additional equilibrium outcomes
x0 = Params.muLearn / ( Params.muLearn + Params.DeltaLearn ) ;                
run MiscFunctions.m

f  = Phi(v,zeta) ;          
w  = NomWagesNoHC(v,zeta) ; % Wages(v,zeta)
wh = NomWagesHC(v,zeta) ;
u  = uRate(v,zeta) ;
H  = MeanHC(x0,Lh,v,zeta) ;
h  = Params.DLearn ./ ( Params.DLearn + Params.phi * u ) ;
s  = Params.delta * zeta ;
ce = Params.ce ;

end

